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Activities in data analysis and numerical simulation of gravitational waves have to date largely 
proceeded independently. In this work we study how waveforms obtained from numerical simulations 
could be effectively used within the data analysis effort to search for gravitational waves from black 
hole binaries. We propose measures to quantify the accuracy of numerical waveforms for the purpose 
of data analysis and study how sensitive the analysis is to errors in the waveforms. We estimate 
that ~ 100 templates (and ~ 10 simulations with different mass ratios) are needed to detect waves 
from non-spinning binary black holes with total masses in the range lOOM© < M < 4OOM0 using 
initial LIGO. Of course, many more simulation runs will be needed to confirm that the correct 
physics is captured in the numerical evolutions. From this perspective, we also discuss sources 
of systematic errors in numerical waveform extraction and provide order of magnitude estimates 
for the computational cost of simulations that could be used to estimate the cost of parameter 
space surveys. Finally, we discuss what information from near-future numerical simulations of 
compact binary systems would be most useful for enhancing the detectability of such events with 
contemporary gravitational wave detectors and emphasize the role of numerical simulations for the 
interpretation of eventual gravitational-wave observations. 



PACS numbers: 

I. INTRODUCTION 

Searches for gravitational waves from coalescing com- 
pact binary systems rely on concrete knowledge of 
the waveform to achieve maximum sensitivity to these 
sources. With LIGO currently acquiring data at de- 
sign sensitivity (see Fig. [T]) , an optimal matched filtering 
search could detect the gravitational waves from binary 
black hole coalescence out to several hundred Mpc. Di- 
rect observation of gravitational waves from these sys- 
tems will have significant and far reaching consequences 
for both gravitational physics and astronomy. 

To date, searches for gravitational waves from com- 
pact binary systems using data taken at LIGO, GEO and 
TAMA observatories have concentrated mostly on binary 
neutron stars and speculative lower mass systems [l[ - 
each element of the binary has a mass below rnj < 3M Q . 
Searches for inspiral waves from higher mass systems such 
as binary black holes and black-hole neutron star pairs 
have used detection templates constructed to match with 
a wide variety of theoretical waveforms 0]. This is the 
first step in searching for one of the most promising and 
tantalizing sources accessible to earth-based gravitational 
wave detectors. As numerical relativity simulations pro- 
duce waveforms, new issues arise when trying to migrate 
this knowledge into the analysis efforts. 

The gravitational waves measured from a compact bi- 
nary system depend on a number of parameters including 



the masses m\ and mi, the spins s\ and s*2, the time of 
merger to, the inclination of the orbital plane t, the phase 
of the orbit at time of merger $0, the distance to the bi- 
nary D, the location of the binary on the sky and the 
polarization angle between the propagation and detec- 
tor frame. For binary black holes these are all the free 
parameters; for neutron stars there are additional param- 
eters that relate to their internal structure and composi- 
tion. In considering the utility of numerical simulations 
in gravitational-wave astronomy, it is critical to under- 
stand the dependence of the numerical waveforms on all 
of these parameters. Some of the dependencies are, in 
principle, in hand already. For example, the functional 
form of the dependence on sky location, polarization an- 
gle and distance is known analytically if the waves can 
be extracted accurately from numerical simulations. The 
time of merger is also easily accounted for by translating 
the computed waveform in time (which is done by apply- 
ing a frequency dependent phase shift to the waveform) . 

Post-Newtonian calculations of the waveforms demon- 
strate the complicated dependence on the other param- 
eters. When all known amplitude and phase terms are 
included in the waveform, it is necessary to explicitly 
measure the masses, spins, inclination and phase. On 
the other hand, the so-called restricted post-Newtonian 
templates (which keep all phase corrections but only the 
leading amplitude term) provide a simplified and efficient 
detection method for non-spinning binaries. The data 
analysis problem (for fixed masses) reduces to the de- 
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tection of a signal of unknown amplitude and phase for 
which the optimal method is well known [3j. 

In this note, we present preliminary investigations 
of waveform accuracy and discuss implications for 
numerical simulations. In Sec. ITTJ we give an overview 
of the form of the gravitational wave strain that is 
measured at detectors, the relationship between this 
and a common way of describing the waveform from 
numerical simulations in terms of Newman-Penrose 
scalars, and summarize the numerical results we will 
use as examples. In Sec. IIII1 we provide a framework 
to obtain an estimate of the accuracy of a particular 
simulation with relation to the data analysis techniques. 
We present examples from the perspective of both 
inspiral and burst searches, and from the inspiral study 
we obtain a crude estimate that, for the purposes of 
detection of non-spinning binaries, a template bank for 
initial LIGO containing approximately 100 waveforms 
will be sufficient; this corresponds to about 10 simula- 
tions with different mass ratios. In Sec. IIV1 we conclude 
with broad discussions of several important and related 
issues, including certain technical aspects of waveform 
extraction from simulations that might introduce spu- 
rious effects — details of this discussion are deferred to 



Appendix [X] We also explain how information from 
current simulations, before fully qualified simulations 
can produce template waveforms, can be used to enhance 
the detection of binary black holes using excess power 
type searches tuned to the right frequencies. Finally, 
we discuss how numerical simulations may be used to 
interpret eventual gravitational-wave observations. Most 
of our analysis focuses on binary black hole systems, 
though in Sec. IIVI we also comment on other compact 
binaries that contain one or two neutron stars. In 
Appendix [B] we provide scaling estimates of the CPU 
time (cost) required to simulate binary systems. This 
should be useful to estimate how feasible various surveys 
of binary black hole merger parameter space would be 
on contemporary computer systems. 



II. GRAVITATIONAL WAVEFORMS FROM 
BLACK HOLE BINARIES 

For a binary black hole system, the gravitational-wave 
strain measured at one of the detectors can be written as 



M*) = 1J ^ PC {^+(0, 5, i, - t .m 1 ,m 2 , si, s 2 , t, $o) + F x (a,5,t,i(;)h x (t - t ; mi, m 2 , si, s 2 , L, $ )} (1) 



where mi and m,2 are the black hole masses, si and s~2 
are the black hole spins, to is the time of merger, the 
inclination of the orbital plane is t, the phase of the orbit 
at time of merger is the distance to the binary is D 
(measured in Mpc) and (a, 5) are the right ascension and 
declination of the binary and ip is the polarization angle 
between the propagation and detector frame. 

Numerical waveforms are usually obtained from the 
Newman-Penrose coefficient ^4. Under a careful choice 
of coordinates, frame and extraction world tube (further 
discussion of which is presented in Scc lIVI and Appendix 
|A"[) the waveform is related to ^4 by 

d 2 

^{h + +lh x ) = ^4. (2) 

It is often convenient to represent the waveform in terms 
of spin- weight —2 spherical harmonics as 

h+, x =J2 cl +,x-2Yi m . (3) 

For concreteness, we consider waveforms extracted from 
numerical simulations performed by Pretorius sim- 
ilar calculations can certainly be done with waveforms 
extracted from other simulations d, 0, B S EH HH • 

The simulations evolve two classes of initial data, a) 
equal mass, quasi-circular co-rotating initial configura- 
tions as calculated by Cook and Pfeiffer [13, HU [H| , and 



b) black hole binary systems formed via the gravitational 
collapse of two boosted scalar field pulses. A detailed 
analysis of the Cook-Pfeiffer (CP) evolutions is presented 
in [l5| . The scalar field collapse binaries (SFCB) have 
non-negligible initial eccentricity, zero initial spin, and 
simulations with mass ratios up to 1.5 : 1 have been per- 
formed (a more detailed description of the equal mass 
scenarios can be found in [1, da]). Note that the scalar 
field is merely a convenient vehicle to create binary con- 
figurations; remnant scalar field energy leaves the vicin- 
ity of the binaries in about one light crossing time (on 
the order of 1/4 orbit), and is dynamically insignificant 
for the subsequent evolution of the binary, and in par- 
ticular the gravitational waves that are generated. Here 
we study three examples from these evolutions — a CP 
d = 16 case (d labels the initial separation between the 
binaries [3), and two SFCBs, one equal mass, the other 
with a mass ratio of 1.5 : 1. The two former evolutions 
exhibit roughly 2.5 orbits before merger, the latter un- 
equal mass case about 1.5 orbits. In all cases the rem- 
nant is a Kerr black hole with spin parameter a ~ 0.7, 
and roughly 3% — 5% of the rest mass energy of the each 
system is radiated away in gravitational waves. Figure [2] 
shows a few samples of the waveforms extracted from 
these simulations, while Fig. [3] demonstrates the conver- 
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FIG. 1: The noise sensitivity curves for the LIGO interferom- 
eters published in June 2006 0] . The blue and red curves are 
the 4km interferometers at Hanford and Livingston, respec- 
tively. The green curve is the 2km Hanford interferometer. 
The LIGO-I noise curve used for the sample calculations in 
this paper is the solid purple line. Compact binaries generate 
gravitational-waves which sweep upward in frequency as they 
inspiral and merge. The frequency (f» 40 Hz) below which 
the noise curve rises sharply determines the longest dynami- 
cal time-scale of the sources to which the LIGO instruments 
are sensitive; this, in turn, translates to a largest mass com- 
pact binary system to which LIGO is sensitive. 

gence behavior of the wave with resolution for the CP 
case. In Fig. [21 data from the highest resolutions avail- 
able are shown. 

Observe that the waveforms depicted in Fig. [2] have 
some noticeable differences. By construction the phases 
all match at t = 0. On axis, the CP and unequal mass 
SFCBs have similar phase and amplitude evolution of 
the waveform several cycles to the left and right of t = 0, 
however, moving toward the orbital plane the similarities 
are less evident. There is also a more rapid decoherence 
between the equal mass SFCB and the other two cases 
moving away from t — 0. One reason for this is that the 
equal mass SFCB example has initial conditions tuned 
to exhibit some "zoom-whirl" type behavior, showing a 
couple of whirl orbits before the system finally merges. 
During the whirl phase the binaries are quite close to- 
gether [inside of what might be considered an inner-most 
stable circular orbit (ISCO)], and moving faster than a 
corresponding point in a quasi-circular inspiral. Hence 
the amplitude of this portion of the wave is quite a bit 
larger than the quasi-circular case, and remains similar 
in magnitude over a couple of wave cycles. These dif- 
ferences are in contrast to quasi-circular equal mass in- 
spiral results obtained by most groups, similar "visual" 
comparisons of which suggest remarkable similarity in 
the waveforms over several wave cycles away from the 
matching point [13, 13- Given that all the latter evo- 



lutions are approximations to essentially the same astro- 
physical scenario, the similarity is not too surprising, but 
nevertheless reassuring. 

It is also interesting to examine the Fourier spectrum 
of these waveforms and to notice the similarities and dif- 
ferences that are manifest. The amplitude of the Fourier 
transform of the gravitational waveform, from the evo- 
lution of Cook-Pfeiffer initial data, shown in Fig. [2] is 
plotted in Fig. |4j This is computed directly from \E r 4(t) 
using the frequency domain equivalent of Eq. ^ which 
gives, for example, 

-4tt 2 / 2 M/) =*?(/) (4) 

where ^>^(f) is the Fourier transform of the real part of 
^4. This avoids introducing artifacts from the numeri- 
cal integration of Eq. To guide the eye, we indicate 
the frequency of the inner-most stable circular orbit es- 
timated by Kidder et al. (l9j 

/ isco ^ 205 ( — J Hz (5) 

and the frequency of qausi-normal ringing given by 

/ qnr « 1600 [1 - 0.63(1.0 - a) 3 ] Hz (6) 

where M is the total mass of the binary, and Mbh arj d a 
are the mass and spin of the final black hole. Notice how 
the power in these waves is predominantly emitted be- 
tween these two frequencies-the initial data is such that 
the binary is orbiting at or near the ISCO frequency. 
Morevover, the spectrum shows a power-law spectrum 
reminiscent of the spectrum of post-Newtonian wave- 
forms. Similar plots are shown in Figs. [5] for the SFCB 
waveforms. It is interesting to note the bump in the spec- 
trum when mi = 7112] this appears to be caused by the 
zooming orbits referred to above. It points to the vari- 
ety of waveform morphologies that might exist for binary 
black hole mergers when eccentricity and even spin of the 
black holes is large. 

The preceding discussion of similarities and differences 
between waveforms is quite heuristic and subjective, and 
thus of arguable merit. One of the main purposes of 
this paper is to propose metrics to both quantify the 
accuracy of simulations and the similarities/differences 
between waveforms, though primarily from the perspec- 
tive of data analysis. The sets of waveforms depicted in 
Figs. [2] and [3l that appear to be significantly different 
due to either numerical resolution effects or different ini- 
tial physics parameters, will provide useful test cases to 
gauge the efficacy of the proposed metrics. 

III. ESTIMATING THE WAVEFORM 
ACCURACY 

The accuracy of numerical solutions is generally deter- 
mined by comparing results obtained at different grid res- 
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FIG. 2: Samples of $4 from the binary black hole merger sim- 
ulations discussed here. Evolutions from three different initial 
conditions are shown: Cook-Pfeiffer d = 16 (CP d=16), and 
two scalar field collapse binaries (SFCB), one equal mass, the 
other with a mass ratio of 1.5 : 1. The top plot shows the real 
part of '3/4 evaluated along the axis 9 — orthogonal to the or- 
bital plane (and azimuthal angle = 0); for brevity we do not 
show the imaginary part as it looks almost identical modulo 
a phase shift. The figures below show the real and imaginary 
parts of ^4 evaluated at 6 — 3tt/8 (note the different vertical 
scale). Here we show both components as there are notice- 
able differences between the two polarizations. In all cases the 
waveform was extracted at a coordinate radius of r = 50m, 
where m is the sum of initial apparent horizon masses; also, 
the time has been shifted so that t = corresponds to the 
peak in wave amplitude, and ^4 has been multiplied by a 
constant complex phase angle to aid comparison. 



olutions. This approach determines the point-wise con- 
vergence of the solutions. In the context of gravitational- 
wave astronomy, only the waveforms themselves are di- 
rectly accessible to observation. It is therefore important, 
if the results of numerical simulations are to be useful to 
gravitational-wave astronomers, to provide a measure of 
the waveform accuracy. 



FIG. 3: A plot demonstrating the dependence on numerical 
resolution of Cook-Pfeiffer d = 16 initial data evolutions. The 
lowest characteristic resolution (dashed line) has a character- 
istic mesh spacing of h, the next lowest one of 3/i/4 (dotted) 
while the finest resolution has a mesh spacing of h/2 (solid). 
The dominant component of the numerical error is in the 
phase evolution of the inspiral portion of the wave. See [l5l ] 
for a detailed discussion of the numerical errors in this set of 
evolutions. 
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FIG. 4: The amplitude of the Fourier transform of the grav- 
itational waveform, from the evolution of Cook-Pfeiffer ini- 
tial data, shown in Fig. [2] The vertical dashed lines are the 
estimated frequency of the inner-most stable circular orbit 
given in Eq. <(Sj and the frequency of the dominant quasi- 
normal mode, assuming a — 0.7, given in Eq. ([6j|. The gray 
shaded region indicates variations in this frequency due to 
10% changes in the mass used. Notice how the power in these 
waves is predominantly emitted between these two frequen- 
cies; the initial data is such that the binary is orbiting at or 
near the ISCO frequency. In addition, the dashed Gray line 
which follows the amplitude is proportional to /~ 5 ^ 6 . While 
this is a convincing fit to the amplitude, we note that there 
is weak evidence for two power laws f~ 7 ^ a , as given by post- 
Newtonian approximations [2$}, below ~ 1.5/jsco and /~ 5 ^ 6 
above that frequency. These simulations do not cover the 
inspiral phase well enough to confirm this result. 
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FIG. 5: The amplitude of the Fourier transform of the gravi- 
tational waveforms shown in Fig. [2] the top panel is an equal 
mass SFCB, and the lower panel is a mass ratio 1.5 : 1 SFCB. 
Note the differences between these spectra and that shown 
in Fig. 2] The bump in the equal mass spectrum arises from 
the hangup of the binary at roughly constant separation the 
a brief whirl phase prior to merger. While these waveforms 
are unlikely to be realized in astrophysical scenarios, these 
spectra provide an example of the rich phenomenology that 
may be realized in binary mergers. 



A. Data analysis formalism 

The standard tool of gravitational-wave data analysis 
is the matched filter. If a signal is present, the detector 
output is 



s(t) = n(t) + h(t;X) 



(7) 



where n(t) is the noise and h(t; A) is the signal. In gen- 
eral, the signal depends on a set of unknown (in ad- 
vance) parameters A. We assume the noise is a zero mean 
(n(t)) = and stationary process, i.e. (n(t)) = and 
(n(ti)n(t2)) — Q(\ti — t 2 \) for some function Q. Here 
(...) denotes the ensemble average over different instan- 
tiations of the noise. 



Define the Fourier transform of the noise by 

/>oo 

n(f) = / n{t)e-' 27r ' lft dt . (8) 



Since the Fourier transform is just a linear transformation 
of the noise time-series, we also have (n(f)) — 0. Further, 
the variance or power spectrum SVid/l) is defined by 

(»(/)»*(/')> = ^n(|/|)*(/-/0- (9) 

Notice that this two point function is diagonal; this is a 
direct consequence of the stationarity of the noise. More- 
over, the power spectrum of the noise will depend non- 
trivially on frequency if the time-domain correlation func- 
tion is not diagonal. Noise with a frequency dependent 
power spectrum is often referred to as colored noise. 

In its simplest form, matched filtering is cross- 
correlating a template waveform w(t; A) with the time 
series s(t) observed by the gravitational-wave detector. 
The matched filter signal to noise ratio (SNR) is 



P(A) 



df 



S(f)w*(f;X) 

Sn(\f\) 



where 



= 2 



df 



Sn(\f\) 



(10) 



(11) 



Notice that the signal-to-noise is normalized such that 
(p 2 ) = 1 if the signal is absent, i.e. h(t; A) = 0. 

Since the signal parameters are not known in advance, 
one must search over all possible values of the parameters 
A to find the template that best matches the signal buried 
in the noise. Then the output of a search (over a small 
chunk of data) will be 



p = maxp(A) . 



(12) 



The largest SNR will be obtained when w(t; A) = 
constant x h(t; A) and A = A. In practice, two important 
issues arise. First, only a discrete set of values of A can be 
searched. This leads to the notion of a bank of templates 
with different parameter values. There is a well devel- 
oped formalism for constructing template banks [2ll | for 
gravitational-wave data analysis. In Sec. IIIID| we will 
comment on the number of templates needed for the bi- 
nary black hole problem and hence the number of accu- 
rate numerical simulations that must be done. Second, 
the theoretical template waveforms may not accurately 
agree with the real signals. It is this issue which we ad- 
dress here. 

The expected SNR is then 



<P(A)> = — 



df 



h(f-\)w*{f-\) 
S{\f\) 



(13) 
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If the template and the signal are the same, then the 
optimal SNR is 



(popt) = — 
Oh 



df 



IM/;A)| 2 



Oh 



(14) 



Define the match [2l| between a waveform and a template 
by 



/' 



(pW) 

(Popt) 



OhO w J_ 00 

The fitting factor [13] 



df 



h(f;X)w*(f;X) 
S(\f\) 



FF = max u, 
x 



(15) 
(16) 

(17) 



is a measure of the distance between a signal and the 
whole template space. If FF = 1, the signal lies inside 
the template space. We will generally use the term match 
to mean something that is maximized over a subset of the 
parameters A and keep fitting factor for the case when all 
template parameters have been maximized over. 



B. Re-parametrization of the numerical templates 

In general, the response of a gravitational-wave detec- 
tor to the waves from a compact binary merger will be 
non-trivial as expressed in Eq. (fT]). Motivated by our 
intuitive expectation that gravitational radiation is dom- 
inated by the quadrupole waves, we have explored the 
following re-parametrization of the numerical waveforms 

, A / 1 Mpc \ r ^ . 

w(t;\) « A I — — — I |cos<Pe + (i - t ;mi,m 2 ) 

+ sin $e x (£ - t ; m 1 ,m 2 )\ (18) 

where A and $ depend on the right ascension a, declina- 
tion S, polarization ip, inclination i and time t. (The vari- 
ation in these constants over the short duration of the sig- 
nal is completely negligible, however the response of the 
instrument to a gravitational wave from a given location 
on the celestial sphere depends on the time of day.) Here 
the plus and cross polarization states e+ iX (t — to', m\ , m 2 ^) 
are just the waveforms extracted on the axis orthogonal 
to the plane of the binary orbit. Here we present evidence 
that this re-parametrization may capture the essential 
features of the merger waveforms sufficiently well for the 
purposes of detecting the waves. On the other hand, the 
full waveforms will probably be needed to extract all the 
possible science. 

To see that this approximation is good enough for 
detection, we re-express the matched filtering SNR in 
terms of these templates. First, notice that the ampli- 
tude A x (1 Mpc/D) cancels out of the SNR defined in 
Eq. (TITJ)) . Hence the SNR can be expressed as 



max 

&,to,mi ,rri2 



z+ cos$ + z x sin<I> 



(19) 



where a w = Da w /(A x 1 Mpc) and 



g(/)g;,x(/;™i.™2) c27ri/t 



(20) 



Sn(\f\) 

For the waveforms considered here, we find that 

.|e x (/;mi,TO 2 )| 2 



M \e + (f; mi ,m 2 )\2 

r MfT) 



df 1 



Sn(\f\) 



(21) 

to better than 3% accuracy. We also find that the two 
polarizations arc almost orthogonal, that is 



df 



e+(f; mi, m 2 )g; (/; m 1 ,m 2 ) 

Sn(\f\) 







(22) 



with typical values ~ 3 x 10 -3 . Hence the normalization 
constant o~ w simplifies considerably and is independent 
of <I> and to: 



df 



|e+(/;mi,TO 2 )| 2 



(23) 



Sn(\f\) 

This allows us to maximize over $ analytically to find 
p= max a~ 1 Jzl(to;m 1 ,m 2 ) + (f ;mi,m 2 ) . 

to ,7711,7712 V 

Below, we will use the plus and cross quadrature matched 
filters to quantify both the accuracy of the numerical 
simulations and the accuracy of the approximation in- 
troduced in Eq. (fig)) . 

Moreover, the match can be written as 



max cos $ + /i x sin $) (25) 

, ,t ,m 1 ,m2 



where 



P+,x 



( z +.,x) 
0~hOe 



(26) 



It is straightforward to see that the validity of the 
re-parametrization in Eq. (| 18[) requires that fi + « 1 
and /j, x ~ 0, independent of the inclination i, when 
h(t,to,mi,m 2 ) = h+{t, to, mi, m 2 , i), and that /x x ~ 1 
and « 0, independent of the inclination i, when 
h(t, to, 77ii, m 2 ) = h x (t, to, mi, m 2 , i). This result is con- 
firmed in Table [U only the mass space needs to be 
searched by the explicit construction of a discrete bank 
of templates. 

It must be stressed, however, that this approach is not 
sufficient for the interpretation of observations and mea- 
surement of parameters. Moreover, there are likely to 
be regimes where this approach is insufficient even for 
detection. 



C. Example of waveform accuracy estimation 

As we mentioned above, differences between the true 
waveform from a binary black hole merger and the tem- 
plate waveform can result in degradation of the SNR. 



TABLE I: The range of matches between waveforms extracted 
at different angles relative to the binary and a template given 
by the waveform extracted at the axis for masses 40Mq < 
M < 200M Q . Note that the ++ and x x entries are all very 
close to unity. The +x entries are close to zero indicating 
that the + and x polarizations are almost orthogonal. This 
suggests that the reparametrization in Eq. (|18p could be good 
enough for detection purposes, with LIGO, in the equal mass 
binary case. 



pp 


h+(to,mi,m2,i) 


h x {to, mi, 1712, i) 


i 


37r/8 


tt/4 


3tt/8 


tt/4 




[0.980,0.995] 
[0.050,0.074] 


[0.990,0.996] 
[0.046,0.069] 


[0.017,0.044] 
[0.989,0.995] 


[0.044,0.066] 
[0.992,0.996] 


SFCB 


h% FCB (t , mi =m 2 ,i) 


/i| FCB (t ,mi =m 2 ,i) 


i 


3tt/8 


tt/4 


3tt/8 


tt/4 


Mx 


[0.981,0.989] 
[0.051,0.075] 


[0.993,0.996] 
[0.024,0.038] 


[0.021,0.044] 
[0.995,0.989] 


[0.019,0.032] 
[0.995,0.998] 


SFCB 


/i + FCB (t ,mi = 1.5m 2 ,i) 


ftx FCB (*o,mi = 1.5m 2 ,i) 


i 


3tt/8 


tt/4 


3tt/8 


tt/4 


A*x 


[0.969,0.973] 
[0.076,0.096] 


[0.986,0.991] 
[0.049,0.061] 


[0.042,0.056] 
[0.987,0.977] 


[0.040,0.050] 
[0.989,0.993] 



In particular, numerically generated waveform templates 
might be inaccurate for any number of reasons, e.g. trun- 
cation errors, instabilities, errors in boundary conditions, 
or incorrect waveform extraction. To get a handle on 
these effects, it is useful to compare waveforms which 
are supposed to represent the same physical process that 
were generated in different ways. 

Here, we present a sample analysis of the waveforms 
presented in |1, [ID, [l6[. Using the formalism outlined 
above, one can compute the matches /i + and /i x for 
the waveform generated at the finest resolution and tem- 
plates at coarser resolutions. This investigation follows 
the standard convergence testing of numerical relativity, 
but with an emphasis on the utility of the waveforms for 
data analysis. Moreover, the answer to this question de- 
pends on the mass of the binary and the detector being 
considered (as it depends on its particular noise curve). 
In Fig. [6l we show the match fi + versus the binary mass 
for the initial LIGO noise curve shown in Fig. [1] For each 
template waveform (i.e. the waveforms from evolutions 
with coarser resolution), there are sets of points: the tri- 
angles indicate the match before maximization over to, 
the circles indicate the match after maximization over 
to. When the match exceeds 0.9, we may be tempted to 
conclude that the finest resolution waveforms are suffi- 
ciently accurate to be used for gravitational-wave data 
analysis. This is not the whole story, however. Referring 
back to Fig. fj] we note that the waves with frequen- 
cies / i5 205(20A/ Q /M) appear to depend on the ini- 
tial data while those above that frequency might reason- 
ably be considered independent of the initial data. Given 
that the current LIGO instruments are sensitive to waves 
above 40 Hz, this suggests that the match should only be 
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FIG. 6: The match between the h+ waveform polarization 
computed as the overlap between the waveforms extracted 
from two different resolution runs. Circles indicate the match 
when the waveforms are allowed to shift in time relative to 
each other. The triangles indicate the match computed as the 
overlap integral between waveforms without the time shift. 
The highest resolution simulation is used as the reference and 
the match with two coarser resolution simulations are com- 
puted. 



trusted for masses M > 100Af Q . 

On the other hand, the small match obtained without 
maximization over time suggests that the evolutions are 
slightly different. In Fig. [3 we show the time shift needed 
to maximize the match. Note that the shift scales linearly 
with mass. This suggests that the difference between 
these waveforms might be captured simply by rescaling 
the mass of the template, i.e. maximizing over both total 
mass and time-delay as one would do in a search. In a 
simple simulation which searched over various template 
masses for a fixed waveform mass, it was found that the 
best match is achieved when the template mass is differ- 
ent from the waveform mass, but the time-shift is similar 
in magnitude to that obtained when using the same mass 
in both the template and the waveform. 

In this concrete example, we have compared the wave- 
forms from the same simulation at different resolutions. 
While the waveforms generated at the finest resolution 
appear to be accurate enough to use as templates in 
searches for gravitational waves, the systematic differ- 
ences between the mass and to which maximizes the 
match hint that the waveforms may not be faithful to 
the physical system. That is, the map between mass of 
the template and mass of the binary (in the standard 
Newtonian sense) may have systematic biases. It will 
be important to explore these issues by comparing wave- 
forms from simulations (starting from the same initial 
data) by different groups. Moreover, detailed exploration 
of the dependence of the waveforms on the initial data 
will also bring information about the faithfulness of the 
wave forms [70]. 
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FIG. 7: The time-offset which maximizes the match (for the 
same mass) between the h+ waveform polarization computed 
as the overlap between the waveforms extracted from two dif- 
ferent resolution runs. The green uses the finest simulation 
as the waveform and the coarsest simulation as the template; 
the red uses the finest simulation as the waveform and the 
next finest simulation as the template. While this shows con- 
vergence, there may still be systematic offsets between these 
waveforms and real physical waveforms. 



D. Comparing waveforms from different 
simulations 

In the previous section, we presented a comparison be- 
tween waveforms extracted from simulations at different 
resolutions in order to gain insight into the accuracy of 
the numerical solutions. Another important step in ex- 
ploring the full parameter space of compact binary in- 
spiral for earth-based detectors is to compare the results 
of different simulations. There are two different reasons 
for making these comparisons: First, comparison of wave- 
forms representing the same physical solution carried out 
using different methods will allow a deeper understanding 
of the numerical issues in this very complicated simula- 
tion problem. Second, it may allow more efficient ex- 
ploration of the parameter space if multiple groups can 
agree on some key test cases and then explore, in detail, 
other regions of parameter space. 

As an example, we compare the waveforms from the 
CP simulations with those from the SFCB simulations 
and between the different SFCB simulations in Fig [5] 
The results quantify the degree to which these waveforms 
are different/similar as seen in Figs. For example, 

the match (maximized over time and phase) between /i^ p 
and the equal-mass SFCB waveforms is greater than 0.95 
for masses M > 200.0. This is consistent with agree- 
ment in both frequency spectrum [/ > 200(40M Q /M)] 
and the waveform. The biggest difference between these 
waveforms therefore appears to come from the eccentric- 
ity of the binary orbit when the black holes first form in 
SFCBs. A similar conclusion holds for the match (max- 



FIG. 8: The match between the h+ waveform polarization 
from one simulation with the waveforms from another sim- 
ulation maximized over both time and phase. The match 

maxt sjn 2 + + H 2 X > 0.95 for masses M > 2OOM0. This gives 
a quantitative measure of the similarities and differences be- 
tween the waveforms presented in Figs. [JHSl 

imized over time and phase) between ft^ FCB (mi = vn^) 
and the SFCB waveform for mass ratio mx/mz — 1-5. 

Finally, we note that the number of templates needed 
to search for gravitational waves from non-spinning bina- 
ries in data from earth-based detectors can be estimated 
as follows. From the CP simulations, we find that the 
match (maximized over time and phase) is fi > 0.97 
for two waveforms with masses differing by « 0.05A/ 
with IOOMq < M < 400M©. Templates along the 
equal mass line would then be laid out with separations 
w 0.1M giving w J*°° dM/{0AM) « 13 templates. If 
we construct a square grid on the mim2-space by draw- 
ing lines through the equal-mass template points, this 
gives !=a 91 templates to cover the square defined by 
50M Q < mi, ?Ti2 < 200M Q . This formula assumes a 
maximum reduction in SNR of 3% due to template mis- 
match. Since a single numerical simulation gives all tem- 
plates along a line of constant mass ratio, this suggests 
that about 10 simulation runs would be needed to cover 
this mass space densely enough for detection purposes. 
This crude estimate should be refined with more detailed 
simulations with different mass ratios and larger initial 
separations. It is important to stress here that this num- 
ber assumes that all the physics and numerical issues 
are under complete control, that no tuning or test runs 
need to be made, and that the waveforms do not get 
much more complicated as the mass ratio changes. All of 
these issues need to be explored before one could gener- 
ate waveform templates that could be confidently used in 
gravitational wave detection. As a result, while the final 
number of templates is not very large, reaching the stage 
where these runs can be made requires much greater ef- 
fort. 
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It is important not to read too much into the esti- 
mate stated above. As we have emphasized, it ignores 
spin and many other important issues in numerical rel- 
ativity. It is also only applicable to earth-based detec- 
tors. The problem is different for higher mass binaries in 
the LISA band. Finally, it only estimates the number of 
simulations needed to enhance detection. As we discuss 
later, this is the first step in using numerical methods 
to extract scientific information from gravitational-wave 
observations, but the larger computational task will be 
extracting accurate information from the data once a de- 
tection is made. 



E. Detectability of numerical templates using 
ringdown niters 

It is illustrative to determine how well existing meth- 
ods of searching for gravitational wave bursts would work 
in detecting numerical waveforms. In particular, we 
would expect that ringdown waveform matched filters 
would work well at detecting the end of the numerical 
waveforms (which do correspond to black hole quasi- 
normal mode ringdown) — especially if it is this portion of 
the waveform that is in the detector band. However, for 
a broad range of masses it is the merger waveform that is 
in LIGO's sensitive band rather than the ringdown wave- 
form: indeed, for the most likely ranges of binary black 
hole masses, the ringdown radiation will be at frequen- 
cies higher than the most sensitive portion of LIGO's 
band. Nevertheless, the numerical waveforms might be 
well-matched by a ringdown template even at frequencies 
below those of the final black hole's quasi-normal mode. 
Will a ringdown matched filter template actually do well 
at detecting these numerical waveforms? Also, will the 
presence of a preceding waveform bias ringdown extrac- 
tion parameters? The answer to both of these questions 
seems to be yes. 

The match, fj,, can be used to measure the ability 
of the ringdown filter to detect the waveform. In the 
case of a ringdown filter the relevant parameters that 
we must maximize in forming the match are the start 
time of the ringdown, to, the central frequency of the 
ringdown, /q, and the ringdown quality factor Q which 
measures the decay time of the ringdown in cycles. A 
ringdown waveform is an exponentially-damped sinusoid: 
cxp(— tt/ot/Q) cos(27t/ t) for t = t — t > 0. As before, 
the match is an indication of the fraction of the signal- 
to-noise ratio that a ringdown filter will obtain compared 
to an optimal filter. The match will depend on the mass 
of the waveform as this determines which portion of the 
waveform is in LIGO's band. From Fig. 0] it can be seen 
that the numerical waveforms do not accurately give the 
gravitational wave energy below 40 Hz for a 100 Mq black 
hole. Since 40 Hz is roughly the low-frequency bound of 
LIGO's sensitive band, care must be taken in interpret- 
ing the match when using low mass waveforms with a 
total mass less than « 100 Mq since a substantial con- 
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FIG. 9: The h+ numerical waveform for a total mass of 50 Mq 
(top), the best-matching ringdown (exponentially-damped si- 
nusoid) filter (middle), and the result of filtering the waveform 
with this filter (bottom), all as functions of time. The match, 
which is the largest absolute value of the bottom plot, is 64%. 
The best match clearly occurs before what we would call the 
ringdown phase of the numerical simulation: this is because 
the ringdown phase is not in LIGO's sensitive band for a total 
mass of 50 Mq — the ringdown filter therefore obtains its best 
match at an earlier (and lower frequency) portion of the nu- 
merical waveform. Although the ringdown filter does not do 
too badly at detecting the numerical waveform, it could not be 
used to measure the quasi-normal mode frequency of the final 
black hole without using information from both the ringdown 
and merger phases. This serves to emphasize the importance 
of developing data analysis techniques which combine infor- 
mation from all three phases of binary evolution [24| . 



tribution to the gravitational wave signal is missing in 
LIGO's sensitive band. 

Figures [9] and [10] illustrate the ringdown filter re- 
sponse to the numerical waveforms for two masses, 
50 Mq and 150 Mq. In these figures there are three 
panels which show (top) the numerical strain waveform; 
(middle) the best-matching ringdown template; and 
(bottom) the result of filtering the waveform with the 
best-matching ringdown template — the match is the 
maximum absolute value of this trace. In both cases 
the match is substantial: 64% for the 50 Mq case and 
70% for the 150 M Q case. In the 150 Mq case (Fig. [TO} 
it is clear that the best-matching ringdown filter is in 
fact matching the ringdown portion of the numerical 
waveform; however, in the 50 Mq case (Fig. [9]) , the best 
matching ringdown waveform is matching the numerical 
waveform considerably earlier than the ringdown phase. 
For the 50 Mq case, the ringdown phase is not in LIGO's 
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FIG. 10: The h+ numerical waveform for a total mass of 
150 Mq (top), the best-matching ringdown (exponentially- 
damped sinusoid) filter (middle), and the result of filtering 
the waveform with this filter (bottom), all as functions of 
time. The match, which is the largest absolute value of the 
bottom plot, is 70%. 



sensitive band. The ringdown filter, while it performed 
reasonably well at detecting the earlier phase of the 
numerical waveform, could not be used to measure the 
quasi-normal mode frequency of the resulting black 
hole. Also note that in the 50 M© case, the match is an 
overestimate because the numerical waveform is missing 
a portion of the late inspiral that would be in LIGO's 
band for this mass, as remarked above. 



IV. DISCUSSION AND CONCLUSIONS 

Recent computational and algorithmic advances in nu- 
merical relativity allow the exploration of strong field 
general relativity in regimes which were hitherto inac- 
cessible. In this paper, we have presented a measure of 
the accuracy of these simulations adapted to the use of 
these results in observational gravitational-wave astron- 
omy. The formalism presented in Sec. MI Al is just the 
match, defined in Ref. [2l[, applied to waveforms from 
various numerical simulations. By applying this metric 
to the equal-mass Cook-Pfeiffer binary black hole simu- 
lations presented in we conclude that these simula- 
tions show convergence within this measure. Neverthe- 
less, we also sound a cautionary note about the match 
as used in Sec. MI Cl to measure the convergence: it only 
measures the convergence, not the physical relevance of 



the ultimate solution. To determine the latter, one must 
also examine a host of other issues including the nature of 
the initial data, the method of waveform extraction, and 
the commonly examined issues of stability, convergence 
and independence on boundary effects. Additionally, the 
currently available waveforms cover just a fraction of the 
relevant sources and analysis similar to those presented 
here will need to be carried out as other cases are treated. 

In the remainder of this section we conclude by dis- 
cussing several outstanding issues in the use numeri- 
cal relativity as a tool for gravitational-wave astronomy, 
including faithful extraction of the waveform from the 
simulations, what could be the most useful information 
that numerical simulations of compact object interac- 
tions could provide in the near-term to enhance the de- 
tectability of gravitational waves, and what information 
from the simulations could help us learn the most about 
compact objects after detection. 



A. Issues in numerical simulations 

Numerical simulations carried out by different codes 
by construction, necessity and available computational 
resources adopt different formulations, employ distinct 
coordinate systems and varied discretization schemes. As 
a result, care is required not only in translating the re- 
sults to data analysis, but also comparing results from 
different codes. Naturally, concentrating on physical ob- 
servables, in our particular case gravitational waves, is 
the sensible way to do this. At the analytical level, a 
well defined approach to compute the radiative proper- 
ties of the spacetime (under natural assumptions) was 
developed in the 60's by taking advantage of future null 
infinity (J^ + )[1S US H3|- However, numerical appli- 
cations dealing with black hole spacetimes can not yet 
reach (this awaits Cauchy-characteristic matching 
or the Conformal equations to be fully implemented in 
the type of systems being discussed here) though there is 
a strong need to compute the radiation produced in the 
problem. To do so within a numerical simulation, two 
approaches are routinely taken. One is based on per- 
turbative methods [HI, 0, H3, HH which relies on a suit- 
able identification of a background spacetime in partic- 
ular coordinates and extracting specific quantities from 
the simulation. A second approach, which has become 
the most common one, makes use of the infrastructure 
developed to calculate the radiation at ^ + but applied 
at a finite distance from the source (see for instance 
[U M, [H M, EE M, M, HI )• While in principle 
this approach can be used beyond the perturbative level 
and is less sensitive to identifying a correct background, 
quantities defined at ^ + need to be translated to finite 
distances where they may not be well or unambiguously 
defined. 

Several key elements, listed below, are in general re- 
quired for faithful extraction of gravitational waves a fi- 
nite distance from the source using the standard result for 
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the relationship between "J 4 and the gravitational wave 
strain (we focus here on items pertaining to J^ + - 
based extraction tools, though similar comments apply 
to perturbative-based methods). For the most part in 
present simulations it is assumed that these conditions 
are satisfied with systematic error less than numerical 
truncation error — of course, these assumptions will even- 
tually need to be verified, and we discuss some sugges- 
tions on how to do this in Appendix [A"l The items listed 
below are not all independent. Furthermore, in theory 
several of the items are not strictly required if during 
wave extraction artifacts induced by "bad" coordinates 
are identified and removed; additional discussion of this 
is also presented in Appendix [Al 

• In the extraction zone the wave travels with unit co- 
ordinate velocity, and the amplitude decays as 1/r. 
If these conditions are not satisfied the extracted 
waves could suffer an error in amplitude and a shift 
in the frequency of the wave (signs of this kind of 
gauge artifact are seen in the CP evolutions with 
generalized harmonic gauge [Hj]). 

• The extraction world tube is assumed to be a ge- 
ometric sphere, and more-over it is assumed that 
the metric of this sphere can be expressed as ds 2 = 
r 2 [d6 2 + sin 2 9d(/) 2 ], where r is the extraction ra- 
dius and (0, 4>) are the usual spherical polar coor- 
dinates mapped onto the extraction sphere. Devi- 
ations from these assumption could, for example, 
lead to artificial mixing of the spherical harmonic 
components of the waveform; and of course the cor- 
rect identification of these harmonic components is 
an important tool in understanding and quantify- 
ing the physics of different merger scenarios. 

• Each point on the extraction world tube is assumed 
to correspond to an inertial observer. Together 
with the above items this is equivalent to assuming 
that the lapse function a = 1 + 0(1 /r) and the shift 
vector induced at the wor ldtubc (3 A = + O(l/r). 
If these conditions are not satisfied all the problems 
mentioned in the preceding items could manifest. 

• Even with a perfect gauge the 0(1/ r) approach to 
Minkowski space could induce spurious effects at fi- 
nite extraction distance; thus the extraction radius 
must be far enough from the source that the 0(1/ r) 
systematic errors in waveforms are smaller than the 
numerical truncation error. This issue can be alle- 
viated in part by a judicial choice of the tetrad used 
to calculate * 4 S& 13, Si ■ 

B. Enhancing the detectability of gravitational 

waves 

Numerical relativity has long been touted as nec- 
essary to doing the best science with ground based 



gravitational-wave detectors since the strongest sources 
of gravitational radiation involve strong gravitational 
fields and the full non-linearity of general relativity. In 
[39j . Flanagan and Hughes laid out the issues relating 
to the detection and measurement of waves from binary 
black holes. They conclude that binary black holes in the 
mass range 25M Q < M < 700M Q may be the strongest 
sources of gravitational waves accessible to earth-based 
detectors. In this mass range, they speculated that most 
of the detectable gravitational-wave energy would come 
from the merger waves emitted between /i SCO and / qnr 
and guessed that about 3% of the binary mass would 
be emitted as gravitational waves from the ringdown of 
the final black hole. Numerical relativity simulations can 
now provide a wealth of information about the merger 
and ringdown phases of compact binaries, even though 
precise connection to the inspiral phase may remain elu- 
sive for some time. 

With the cautionary note that so far numerical sim- 
ulations have provided input on a (small) subset of the 
physical parameter space and that, in particular, spin- 
orbit interactions might strongly influence the modeled 
waveforms, valuable insights can be drawn with the cur- 
rent knowledge. Consider the simulations of CP initial 
data discussed here. We can immediately make sev- 
eral qualitative observations about the waves (see also 
the relevant discussions in [IB]). First, the waves sweep 
smoothly upward in frequency from /i sco to / qnr ; dur- 
ing this phase the time domain amplitude also increases 
monotonically. About 3-5% of the binary mass is radi- 
ated in the last plunge-orbit, merger and ringdown. The 
analysis in Sec. IIIIEI suggest that upwards of ~ 70% 
of the wave energy can be attributed to the ringdown 
waves (though some portion of the wave matched by 
the ringdown template could be associated with the late- 
inspiral part of the wave). Perhaps most important for 
detection of the waves is that the bulk of the wave en- 
ergy is emitted in a time-frequency volume defined by 
102.5(40M Q /M) < / < 5OO(4OM /M) within a time 
duration At « 0.03(M/40M o ). In the language of An- 
derson et al. |4(j, the time-frequency volume of the sig- 
nals is V — 12 with a definite mass-scaling for the fre- 
quency band. Figure 1111 compares the sensitivity of a 
matched filtering search to an excess power search for 
several different time-frequency volumes. In the context 
of these equal mass, non-spinning binaries, it suggests 
that a naive excess-power search would miss about half 
of the signals detectable by a matched filtering search for 
the merger only. Nevertheless, this information, when 
carefully combined with a search for gravitational waves 
from the inspiral and ringdown phases of binary evolution 
could provide a near optimal search for these sources. 

Dynamical simulations of binary black holes have so 
far only been carried out for a very limited number of 
parameters (e.g. mass ratio and black hole spin), and 
it will be very important to systematically investigate 
the dependence of the waves on these parameters (com- 
pare [3, Sl|). For example, black hole spin may effect 
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FIG. 11: The relative effectiveness r\ of an excess-power 
search with known time-frequency volume V compared to 
a matched filtering search with N effectively independent 
templates assuming both searches are tuned to give a 1% 
false-alarm probability. The relative effectiveness is defined 
as r\ = Dgg%/D^% where -Dgg% is the effective distance at 
which 99% of the sources are detectable in each search. Given 
the estimated time-frequency volume of the merger signal 
from equal mass binaries, i.e. V ~ 12, this suggests that 
a sub-optimal search strategy could still detect about half the 
sources that could be detected in a matched filtering search; 
simple enhancements to the excess power method should be 
able to do better. This data is taken from Fig. 4 of Ref. [40l |. 



it will be impossible to populate a template bank solely 
with the results of numerical simulations in the near to 
medium term. It is therefore interesting to devise approx- 
imate methods which might capture the essential features 
of the merger phase. One possible approach would be to 
use PN methods for the early phase and close limit ap- 
proximations for the late phase with a judiciously chosen 
behavior in between. Adopting such approximate meth- 
ods requires a careful understanding of the sensitivity 
of detection methods to differences in the approximate 
waveforms. 



C. Learning about compact objects through 
gravitational-wave observations 

While the direct observation of gravitational waves will 
be a huge achievement, we hope that the first detection 
will only mark the beginning of gravitational-wave as- 
tronomy. This field is built on the premise that we can de- 
code information about the sources of gravitational waves 
from the signals observed at a detector. To achieve this 
goal, we need the ability to simulate the generation of 
gravitational waves by various sources. Moreover, it is 
the imprint of the sources on the waves that will carry 
some of the most interesting information. 



1. Binary Black Holes 



both the amplitude and the duration of the merger waves. 
For aligned black hole spins the merger takes longer than 
for anti-aligned black hole spins [4JJ. In addition, the 
interaction between spin and orbital angular momentum 
may play a role in the dynamics of the merger (even 
though the speculations of [12] have not been confirmed 
by the dynamical simulations of [Hj]). 

An important next step in confirming these specula- 
tions is to use the waveforms from numerical simula- 
tions as sample signals in real detector noise and passed 
through the current detection pipelines. This would be 
facilitated by the development of an archive of gravita- 
tional waveforms in a uniform format that could be used 
openly by the gravitational-wave detection community to 
calibrate their searches and their pipelines. Activities in 
this direction are already under way [44| ; with the inten- 
tion of developing into a useful resource for gravitational- 
wave astronomy. As an aside, to expedite their use in 
data analysis, these waveforms should be provided as 
equally sampled time series, with time in units of sec- 
onds and scaled to a physical distance of 1 Mpc from the 
binary. 

Despite recent computational and algorithmic ad- 
vances, numerical simulations are costly and will be for 
years to come. In Appendix [Bl we present order of mag- 
nitude estimates of the computational cost of simulations 
needed to produce parameter space surveys of a given ac- 
curacy and physical evolution time. This suggests that 



By exploring the results of merger simulations from 
the perspective of data analysis we can obtain a better 
picture of what we would be able to learn about compact 
object interactions from future observations. Therefore 
it will be very important in the near future to perform 
surveys of a wide variety of initial data parameters (in 
particular varying mass ratios and spin vectors), not so 
much to build template libraries, but to understand what 
the broad features of merger look like through the lense 
of a gravitational wave detector. For example, it has long 
been anticipated that the onset of the binary merger, at 
which the slow and adiabatic binary inspiral changes into 
a dynamical plunge, would occur at an inner-most stable 
circular orbit and would leave a characteristic signature 
in the gravitational wave signal, which could be measured 
in gravitational- wave detectors. While the recent dynam- 
ical simulations of binary black holes do not reveal any 
abrupt change in the waveform, there does appear to be 
a break in the frequency spectrum of the waves which oc- 
curs somewhere between the predicted ISCO frequency 
and the quasi-normal mode frequency — see Fig. [H and 
further discussion of this in [l5| . It is certainly plausible 
that this break in frequency becomes more pronounced 
as one moves away from the equal mass, non-spinning 
regime. 

Black hole coalescence has also been regarded as giving 
rise to an arena where the strong-field, non-linear regime 
of general relativity will be clearly revealed to observers. 
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The early simulations show perhaps a disappointing lack 
of such features, where except for a very short and 
smooth transition between inspiral and ringdown, much 
of the waveform can be understood using perturbative 
techniques. Another way of stating this is that all wave- 
forms to date are dominated by the quadrupole har- 
monic. In general one would expect non-linear effects 
to result in mode coupling. Again, that we do not see 
significant higher order harmonics could be due to the 
restricted initial conditions so far considered; however, 
at the very least this is telling us that manifest strong- 
field effects are not ubiquitous in this type of black hole 
collisions, and the community will need to search harder 
to find richer regions of astrophysically relevant param- 
eter space. Whether this observation remains as such in 
more generic cases will have strong consequences for the 
simulation and analysis sides. On one hand, decipher- 
ing the non-linear effects would require significantly more 
accurate simulations and a considerably denser template 
bank for data analysis. On the other hand however, these 
templates could be parametrized in a rather simple form 
like that in Eq. ([15]). 



2. Binary Neutron Stars 

For binary neutron stars the situation is significantly 
different from that for binary black holes. Depending 
on the equation of state, stellar masses and spins, the 
merger of binary neutron stars may be triggered either 
by a plunge after the two stars reach an ISCO, or by 
Roche lobe overflow (see, e.g. [iHlifjj]). In either case, the 
merger is expected to occur at a gravitational wave fre- 
quency of approximately 2 kHz, outside of LIGO's most 
sensitive regime. That means that the current LIGO con- 
figuration is more sensitive to the inspiral phase, which 
may be well approximated by post-Newtonian calcula- 
tions, than the merger of binary neutron stars, which 
has to be modeled with numerical relativity. Hence, the 
role of numerical relativity in observing binary neutron 
stars is different from its role for binary black holes. In 
the near-term, numerical simulations could be used to 
validate the post-Newtonian approximation to the wave- 
forms in the LIGO frequency band; whether this is possi- 
ble on current generation computing facilities remains to 
be seen. Furthermore, numerical relativity may provide 
guidance for the design of future configurations, given the 
astrophysical scenarios that seem particularly promising. 
In the long-term, the prospect of observing gravitational 
radiation from the merger of binary neutron stars is very 
exciting because it is very rich in physical effects that may 
play an important role. Unlike binary black holes, which 
are governed entirely by Einstein's field equations, the 
dynamical evolution of binary neutron stars also depends 
on the equation of state, magnetic fields, radiation and 
neutrino transport, and possibly other effects. Realistic, 
nuclear equations of state have already been adopted in 
simulations of binary neutron star mergers pTj , and nu- 



merical codes that incorporate general relativistic mag- 
netohydrodynamics have been developed (e.g. [H, [i^]). 
Detecting a binary neutron star merger may therefore es- 
tablish important observational constraints on these as- 
pects, in particular the equation of state. Clearly, this is 
a very exciting prospect. 

As discussed above, it is unlikely that the current 
gravitational-wave detectors could observe the details 
of the merger. Numerical relativity may nevertheless 
play an important role for the purposes of data analy- 
sis, namely by identifying features in the gravitational 
wave signal that could provide particularly important in- 
formation. A concrete example is the question whether 
or not the merger remnant promptly collapses to a black 
hole. Depending on the equation of state, binary neutron 
stars may either collapse to a black hole on a dynam- 
ical timescale after merger or the remnant may form a 
"hypermassive" neutron star supported against collapse 
by virtue of differential rotation [47l l50f . For binaries 
with fixed masses the pre-merger gravitational- wave sig- 
nal is quite similar, but the post-merger signal differs 
significantly for the two scenarios. Distinguishing these 
post-merger signals would therefore provide an important 
constraint on the neutron-star equation of state. Unfor- 
tunately, the post-merger signal typically has a frequency 
close to 4 kHz [47j ] and may require advanced, perhaps 
specialized, detectors to extract significant information 
(see also the discussion in 5l|). Investigations of these 
waves and the information they carry could therefore be 
used to inform development of advanced gravitational- 
wave detectors. 



3. Black hole-neutron star binaries 

For mixed black hole-neutron star binaries the situa- 
tion is different again. The inspiral can have two very 
different outcomes: either the neutron star is tidally dis- 
rupted by the black hole before it reaches the ISCO, or 
it reaches the ISCO first and spirals into the black hole 
more or less intact. A very crude scaling argument sug- 
gests that the separation of tidal disruption stid is ap- 
proximately 



gtid 
Mbh 



2/3 



-Rns 



/MnsY 
VA/bh/ M NS ' 



while the ISCO is located approximately at 



sisco 



Mi 



G. 



BH 



(27) 



(28) 



If Stid > siscOj the neutron star is disrupted before it 
reaches the ISCO, and vice-versa. Since most neutron 
stars are expected to have a mass of slightly more than a 
solar mass and ratio -Rns/M^js ~ 5, the outcome mostly 
depends on the black hole mass Mbh (see Ref. (52j|). 

The neutron star can only be disrupted tidally out- 
side the ISCO for relatively small stellar-mass black holes 
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with Mbh ;$ 5AfNS- This regime is very interesting for 
a number of astrophysical reasons. Such a disruption 
may act as the central engine of short gamma ray bursts, 
and, as for binary neutron stars, a detection may pro- 
vide useful constraints on the equation of state. It is also 
this regime that requires numerical relativity simulations 
for quantitative predictions. Extending the above crude 
estimates, the orbital frequency at the onset of tidal dis- 
ruption is approximately 



Mi 



BH 



3 tid 



1/2 



M- 



NS 



R 3 



1/2 



(29) 



which is the same order as the inverse of the neutron 
star's dynamical timescale. As for the merger of binary 
neutron stars, the corresponding gravitational-wave fre- 
quency is outside of LIGO's most sensitive wave band. 
The role of numerical relativity is therefore again to ex- 
plore these scenarios and identify particularly interesting 
frequency regimes. 

Numerical relativity simulations of mixed black hole- 
neutron star binaries are not as far advanced as those 
of black hole binaries or neutron star binaries. So far, 
fully relativistic, self-consistent studies exist for quasi- 
equilibrium models [53], HU, (see also [56] ) . The first 
fully self-consistent, dynamical simulations of the binary 
inspiral and the tidal disruption of the neutron star have 
been announced very recently [53] (see also [6(| for simu- 
lations of extreme mass-ratio binaries within the so-called 
Wilson-Mathews approximation, which assumes that the 
spatial metric remains conformally flat). Othe r gr oups 
have also initiated studies of mixed binaries [5a. |59j. 
Clearly, more comprehensive studies of tidal break-up in 
black hole-neutron star binaries remain an important and 
urgent goal. 

These calculations will help to address several very im- 
portant questions. One such question is whether the tidal 
disruption of a neutron star in a mixed binary may lead to 
an accretion disk that is large enough to power a gamma 
ray burst. Another question concerns the nature of the 
mass transfer, which could proceed on a dynamical or sec- 
ular timescale, or could even be episodal (see, e.g., the 
discussion in [601 ] and references therein). The nature of 
the tidal disruption has great impact on the gravitational 
wave signal emitted in the process. Presumably the dis- 
ruption will leave a signature in the signal at a frequency 
related to orbital frequency. As inferred from equation 
(|29p . this frequency carries important information about 
the neutron star and its structure. A quantitative under- 
standing of these issues clearly requires detailed numeri- 
cal relativity simulations. 
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APPENDIX A: CALCULATION OF RADIATION 
AND SYSTEMATIC EFFECTS 

In this section we discuss in more detail possible sys- 
tematic effects resulting in the calculation of waveforms. 
We concentrate in particular in the approach based on 
the Newman-Penrose (spin-weighted) scalar ^4 as is the 
most commonly employed, however similar issues arise in 
the perturbative approach as well. 



1. From the analytical to the numerical arena. 

^4 is a particular combination of the Weyl tensor in a 
suitable frame and coordinates. Its leading order ^4 in 
a suitable Taylor expansion off future null infinity {J? + ) 
provides the gravitational waves. From now on, as we 
concentrate on the extraction of gravitational wave we 
will drop the supra-index "°" from all related quanti- 
ties, though it must be understood that we are referring 
to the leading behavior. A related quantity, the shear 
of the outgoing null rays a plays a crucial role in the 
calculation of radiated momentum and angular momen- 
tum. A few key features are trivially satisfied at J? + 
by construction which makes the unambiguous calcula- 
tion of the radiative properties of the spacetime possible. 
Namely: the metric is exactly flat, the location of the ex- 
traction worldtube is unique (up to time u translations), 
and a powerful structure (like the asymptotic transfor- 
mation group, asymptotic flatness condition and peeling 
behavior) allows for a clean definition of the sought-after 
quantities. This ensures a unique calculation of the ra- 
diative properties of the spacetime. 

Unfortunately, most codes are unable to calculate these 
quantities at future null infinity and thus $4 is obtained 
at finite distances fn\. Since the extraction worldtube 
does not represent a flat surface, key ingredients are miss- 
ing which can have a non-trivial impact in the calculated 
quantities, even when the suitable decay of this quantity 
is materialized (i.e. peeling is satisfied). For instance, 
the suitable frame allowing for the calculation of radia- 
tion at future null infinity - known as a Bondi frame - is 
such that the angular part of the metric is exactly that of 
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the unit sphere (i.e. the angular metric induced at any 
given time is exactly that of the unit sphere), there is 
no induced 'shift' in the coordinates of observers at 
{9uA — 0) and the observers retarded-time u is affincly 
parametrized (g ur = 0). These conditions play a cru- 
cial role in several aspects: (i) inertial observers maintain 
constant angles and clocks that tick at a constant rate (ii) 
the variation of the angular part of the metric near ^ + 
is solely due to gravitational waves and (iii) the simple 
relation $4 = a uu holds (with a denoting the complex 
conjugate of a). This relation, together with the Bianchi 
identities at future null infinity, is employed to replace 
the appearance of a by suitable integrals of ^4. 

At the numerical level, the worldtube is routinely de- 
fined by a Cartesian timelike worldtube at x 2 = r 2 inter- 
secting hypersurfaces at t = const. The induced metric 
on this worldtube generically does not satisfy the con- 
ditions qab — t 2 1ab + 0(r) (with qAB the unit sphere 
metric, <?tt_= 1 gtA = 0. Consequently, as discussed in 
detail in [6j|, if 'J 4 ^ a (with a " ' " indicating <9 t ), 
observer's clock rates (at different locations on the ex- 
traction worldtube) tick at dissimilar rates and do not 
stay at constant angular locations. These issues can in- 
troduce systematic effects which can affect the predicted 
outcome. For instance, 

• If ^4 ^ fj, commonly employed formulae which re- 
places er° by integrals of ^4, lack non-trivial con- 
tributions from products of a and (time deriva- 
tives of) the conformal factor F relating the an- 
gular metric gAB to the unit sphere metric qAB ~ 
see Eq. (|A1|) . Recall that since gAB is the metric 
of a sphere, it is conformally related to that of the 
unit sphere by gAB = r 2 F 2 q A B 

• If gtA 7^ 0, inertial observers suffer a rotation (from 
the induced shift Pa — gtA)- Therefore, predictions 
like the waveforms at a particular angle will be af- 
fected as inertial coordinates are shifting around 
the worldtube. This will influence the extraction 
of multipole contributions, since the spin-weighted 
spherical harmonics employed in such a task do not 
take into account the shift in the angular coordi- 
nates (see [gil for a related discussion of these is- 
sues at 

• If gtt 7^ 1 the radiation measured by different ob- 
servers at a constant t — const slice does not corre- 
spond to the same inertial time. Thus the extracted 
waveforms would have to be mapped, at each angle 
to the real inertial time. 

It is clear that these issues can introduce systematic 
effects that can be either corrected or at least estimated 
within a given simulation. Notice that these issues can 
not be completely addressed with an improved tetrad 
choice as any tetrad must satisfy g a b = 2l(a n b) ~ ^^{a m b) ■ 
Thus, the induced metric in all cases will be the same. 
Nevertheless a convenient choice of tetrad aids, in par- 



ticular, to alleviate issues related to the proximity of the 
extraction worldtube to the source. 

It is useful to first at least determine the magnitude 
of the effects the issues discussed above might have in a 
given simulation. Then if required, these effects can be 
corrected by suitably modifying the employed expressions 
to remove many of the ambiguities. In what follows we 
describe the first part, while a discussion of the second 
will be presented elsewhere [63j and applied in relevant 
scenarios [65], 

First and foremost is estimating the amount that the 
above mentioned effects can have on ^4. We will ignore 
here those coming due to the extraction at finite distances 
as a comparison of obtained results at different radii can 
be employed to estimate this effect. Furthermore, assum- 
ing the extraction takes place sufficiently far away, the 
difference between a derivative in the time labeling the 
timclikc hypersurfaces and null hypersurfaces -beyond 
gauge factors which we will discuss later- is given by 
the contribution of spatial derivatives off the extraction 
worldtube times 1/r factors due to an approximate "po- 
tential" from the source. This effect is also controlled by 
placing the extraction sphere sufficiently far and compar- 
ing the obtained results at different radii. Under these 
assumptions, the main contribution to the error in as- 
signing the radiation to the straightforwardly calculated 
^4 is given by the discrepancy of the extraction sphere 
being in agreement with that of a Bondi frame. This 
is made evident by the failure of the induced metric on 
the sphere to be that of the unit sphere (66|. In turn one 
has gAB = t 2 SaBi which can be re-expressed in terms of 
Sab = F(t, 0, <p) 2 qAB with qAB the unit sphere metric. 
The correction to ^4 is given by 



*2 = & 



*^M-*m (ad 



with 9 the eth operator which is a particular combination 
of derivatives on the sphere [22j. Thus, unless F = 0, a 
non-trivial correction must be considered from the fact 
that the extraction sphere is not inertial. Certainly this 
will occur as the very radiation one is trying to com- 
pute will be responsible for accelerating the sphere. The 
key message is then to estimate the role F will play. 
To this end two calculations can be performed. First, 
a partial answer on the value of F can be obtained by 
simply evaluating F = det(g ab) / 'det(r 2 'qAB) ■ If F 7^ 1, 
then F ^ 1 and it can therefore play a non-trivial role. 
However, if F = 1 it is not necessarily the case that 
F = 1. A more involved, though now complete, recipe 
to obtain F can be easily obtained by computing the 
Ricci scalar associated with gAB and qAB and the fact 
that the two metrics are conformally related by F. This 
gives rise to the expression 1Z = 2 (F 2 + V^V^ logF). 
Notice that F = 1 is a solution if gAB is the unit 
sphere metric, hence, short of obtaining a solution for 
F, an estimate of ignoring this fact can be obtained 
by Ef = \\~R- — 2||. Naturally, if Ep remains well be- 
low the measured waveforms, the straightforward use of 
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"J 4 would be warranted. To simplify the numerical cal- 
culation of 1Z one can make use of the Gauss-Codacci 
relations for a 2-dimensional hypersurface S in a three- 
dimensional manifold £ and employ quantities readily 
available on X. Defining the extrinsic curvature of S by 
n a b = h c a hfS/ c Sd (with h a b = j a b-s a Sb, the induced met- 
ric on S with normal s a — V a i?, s a — s a /(s a s a ) 1 / 2 and 
V c 7ab = 0), a straightforward calculation indicates 



K = (3) i?-2 (3) i? a , 



a „b 



bS S 



cd 



(A2) 



In addition to the conformal factor F being taken prop- 
erly into account, a Bondi frame satisfies that observers 
measure an affine time along ^ + and proceed along con- 
stant angular coordinates. These conditions will be met 
unless gtt = l,9tA = 0. Here again norms could be de- 
fined as Ecu = \ \{gtt - 1)||, &GtA = \ \gtA\\ so as to obtain 
an estimate of the effect these issues might have in the 
extraction process. 



2. Coordinate conditions and extracted 
quantities — an example 

To illustrate some of the effects of coordinate condi- 
tions that are not well adapted to the extraction mecha- 
nism we adopt a spacetime containing linearized gravita- 
tional waves (67[. For our particular example, the space- 
time is described in terms of the following line element 

ds 2 = -dt 2 + (1 + Af rr )dr 2 

+2Bf rg drd6 + 2Bf r4> sm(9)drd(j) 

+ ^ + Cfle + Af 2 ee )r 2 de 2 

+2(A - 2C)f e<t> r 2 sm(6) 2 d9d<l> 

+ {l + Cf^ + Af 2 ^r 2 S m(e) 2 d^ 2 (A3) 

where 

f rr = sin(6>) 2 cos(20) ; f r g = sin(20) cos(2</>)/2 
f r< l> = - sin(6>) sin(20) ; f e<j) = cos(6») sin(2</>) 
fee = -fU = (l + cos(0) 2 )cos(20) 
/| fl = -cos(20) ; f H = -co S (6) 2 cos^) 



and 



.4 



B 



-3sin(Y") 9cos(y) 9sin(F) 



cos(F) + 3 sin(Y) 6 cos(Y) 
6sin(r)' 



(A4) 
(A5) 
(A6) 



c = 1 ^ sm(y) _ 2cos(r) _ 9sin(r) (A?) 



21 cos(F) 21 sin(F) 



with Y = t — r and for simplicity we have adopted F — 
s'm(t — r) in reference [67|. While a detailed discussion 
of of problematic issues that can arise in the extraction 
proceess as well as ways to handle them will be discussed 
elsewhere [63[ we here illustrate the effect that some of 
these will have in a simple scenario. 

Notice that the line element IA3I satisfies all the men- 
tioned properties, in particular g u — 1, gtA — and 
9AB = r 2 q A B + 0(r) (with A = B,4> and q AB = 
diag[l, sm(9) 2 ]). For clarity, we will concentrate on two 
very simple cases that will violate these conditions de- 
fined by r — > g(t)r and <f) — * (p+Lut.The former introduces 
a time-dependent variation in the location of the extrac- 
tion radius with respect to a physically defined areal ra- 
dius. The latter induces a shift along the direction. 
Following the commonly used approach one obtains for 

^ = ^mp^n cos{2 ^ +LOt)) (A9 ) 



4ff 

c^oj = zS in(Y)cos W „ ; 



2.9 



sm(2{(/) + ujt)) (A10) 



with Y = t — rg; clearly, while a calculation of the ra- 
diated energy would be immune to the value of u> the 
actual waveforms would be affected. Additionally, both 
would be affected by the functional dependence of g(t), 
which even when chosen initially to be unity, it will vary 
as waves propagate outwards affecting the spacetime. 
Clearly a generic situation would not be as simple as the 
one analyzed above, though it serves to make evident how 
these issues can be obscured in a non-trivial extraction 
of physical quantities in an otherwise correctly obtained 
numerical solution. 



APPENDIX B: ESTIMATING THE 
COMPUTATIONAL COST OF BINARY BLACK 
HOLE MERGER SIMULATIONS 

Here we present an order of magnitude estimate of 
the computational cost of simulating a binary black hole 
merger for a specified number of orbits, and where we 
want to extract some feature of the solution to within 
a given accuracy. More precisely, if the CPU run time 
and estimated net solution error of a fiducial simulation 
are known, we give an order of magnitude estimate for 
the CPU run time that a second simulation will take to 
obtain a new result with the same accuracy but from 
an evolution that completes a different number of orbits 
before merger. 

Assume the numerical error E(t, h) in some desired 
property of the solution (for example the phase evolution 
of the gravitational waveform) has the following depen- 
dence on physical time t and characteristic discretization 
scale h: 



(A8) 



E(t, h) oc t q h r ' 



(Bl) 
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Here q is a positive constant of order unity, and m is the 
order of the discretization method. In general the growth 
of error will be more complicated than this simple power 
law, though for an order of magnitude estimate this ex- 
pression is sufficient. Note also that we have assumed 
the code has been "cured" of any exponential growth in 
error. In an adaptive code there will be several mesh 
spacings h, though the scaling relationships derived be- 
low will still be correct if all mesh resolutions are changed 
by the same factor when the resolution is changed. Equa- 
tion (|B1[) is strictly only valid for finite-difference codes, 
though in the final expression below we can take the limit 
as m — > oo to get estimates for spectral codes. This will 
give us an idea of the scaling of present pseudo-spectral 
codes in the regime where the leading source of error 
comes from the spatial discretization and not the finite- 
difference time stepper. 

A more useful parameter describing the physical run 
time t of the simulation is the number of orbits n com- 
pleted before coalescence, and we will assume that the 
inspiral regime of the merger dominates the run-time. 
We can use the leading order Post Newtonian expression 
for equal mass, quasi circular inspirals to estimate n(t): 

n(t) oc t 5 / 8 (B2) 

The final ingredient in our scaling estimates will be 
the manner in which computational run time T scales 
with mesh spacing h for any optimal grid-based solution 
method of a 3 + 1 dimensional system of partial differen- 
tial equations 

T(M)«^ (B3) 

The first question we can now answer is the follow- 
ing: given that simulation A required Ta CPU hours to 
complete a simulation of a binary system exhibiting ua 
orbits before coalescence, and from which we extracted a 
desired quantity with error Ea, how long Tb will it take 
to run second simulation of a similar binary system, now 
with n b orbits before coalescence but the same net error 
E B = E A 7 Using (|B1IIB3I) we find 

/„ \ (8/5+32g/5m) 

T b =Ta(^A . (B4) 

For example, assuming linear growth of error (q = 1), 
we get T B = T%, where z = 4.8, 3.2, 2.7, 2.4, 1.6 for 
to — 2, 4, 6, 8, oo respectively. Note that the differ- 
ence going from 2 nd to 4 th order accuracy is quite signif- 
icant, as is the jump from 4 th to "spectral" convergence 



to = oo. The same holds when estimating the accuracy 
achieved by different order-of-accuracy operators in mod- 
eling modes in the solution (68l . [69^ . This implies that 
while higher order methods are important, a significant 
gain is already achieved at 4 th order. 

An application of the preceding expression is to the 
recent survey of unequal mass inspirals presented in [HI ] . 
The 4 th order accurate code (to = 4) discussed there 
is quite fast by today's standards, and they were able 
to perform the survey utilizing a total of about 150, 000 
CPU hours (recall however that the cost of these simu- 
lations is in practice further alleviated by exploiting the 
problem's symmetry). The majority of initial conditions 
ran exhibited about 2 orbits before merger. Suppose the 
survey were repeated (including calibration runs, etc.), 
but now starting with initial conditions resulting in 4 or- 
bits prior to merger. Equation (|B4[) suggests (assuming 
q = 1) it would take around 1.4 million CPU hours to 
complete with the same level of overall accurac y. E arly 
comparisons of numerical versus PN waveforms [15| sug- 
gest more than 4 orbits are needed to begin to study 
the adequacy of various PN approximants. Suppose 10 
were sufficient, and the survey of [1(| were repeated for 
the purpose of PN comparisons — (|B4[) then suggests 170 
million CPU hours would be needed. This would require 
a 20, 000 node cluster running continuously for 1 year. 
Of course, we are not suggesting that such a survey is 
necessary for PN comparison purposes, we are merely 
illustrating (|B4p using actual data for the reference sim- 
ulation "A" as given in 10]. 

A second interesting question we can give a rough an- 
swer to is, given simulation A has accumulated an error 
Ea for an Ua orbit simulation, what is the expected er- 
ror for simulation B that completes ub orbits, assuming 
simulations A and B have identical resolution Ha = ^s?: 

fn B \ {8q/5) 
E B = E A . (B5) 

At a first glance it might seem strange that the accu- 
mulation of error is independent of the order of the dis- 
cretization scheme. Though observe that this does not 
imply that a lower order scheme is just as "fast" as a 
high order method even if both methods exhibit similar 
growth factors q, for in general it will take more resolu- 
tion for a low order method to get to the same level of 
error E as a higher order method. From (|B3|) the ratio 
of run times Th t /Th 2 for two different resolutions hi and 
/i2j regardless of the order of convergence for an optimal 
solution method, scales as (ft.i//i2) 4 - 
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